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ABSTRACT 

An ultrarelativistic relativistic study of axisymmetric Bondi-Hoyle accretion onto a 
moving Kerr black hole is presented. The equations of general relativistic hydrody- 
namics are solved using high resolution shock capturing methods. In this treatment 
we consider the ultrarelativistic limit wherein one may neglect the baryon rest mass 
density. This approximation is valid in the regime where the internal energy of the 
system dominates over the rest mass energy contribution from the baryons. The pa- 
rameters of interest in this study are the adiabatic constant T, and the asymptotic 
speed of the fluid, Voq. We perform our simulations in three different regimes, subsonic, 
marginally supersonic, and supersonic, but the primary focus of this study is the pa- 
rameter regime in which the flow is supersonic, that is when Voo ^ cf . As expected 
from previous studies the supersonic regimes reveal interesting dynamics, but even 
more interesting is the presence of a bow shock in marginally supersonic systems. A 
range of parameter values were investigated to attempt to capture possible deviations 
from steady state solutions, none were found. To show the steady state behaviour of 
each of the flows studied we calculate the energy accretion rates on the Schwarzschild 
radius. Additionally, we also find that the accretion flows are dependent on the loca- 
tion of the computational boundary, that if the computational boundary is located 
too close to the black hole the calculated flow profiles are marred with numerical ar- 
tifacts. This is a problem not found in previous relativistic models for ultrarelativistic 
hydrodynamic systems. 

Key words: accretion, accretion discs - black hole physics - hydrodynamics - rela- 
tivity - shock waves - methods: numerical 



1 INTRODUCTION 

During the radiation dominated era the universe is thought 
to be comprised of radiation fluid, that is to say a fluid 
with internal energy which is sufficiently large that the 



rest mass density is negligible (Miller 2005 Musco et al 



2005 Musco 2008 1. During this phase primordial black holes 



(PBH) would be created ( Carr and Hawking 1974 1 and grow 
due to the accretion of this radiation fluid Sivaram 2001 1. 
Over the last several decades there has been interest in the 
formation of primordial black holes (PBH). In particular 
there has been focus on the question about the growth rate 
of the PBH. Recent studies of primordial black holes by |Cus-| 
todio and Horvath (2002) suggest that their accretion rates 



are relatively low, so that the mass of the black holes does 
not increase much through the accretion process. 

Previous studies used a simplified Bondi-type model 
(stationary accretor) with an ideal gas equation of state 
to show that the accretion rates would not be sufficient 
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to appreciably modify the mass of the PBH if the hori- 
zon of the PBH is smaller than the cosmological horizon 



(Zel'dovitch and Novikov 1967 Carr et al. 20101. In the 



event that the PBH is on the order of the cosmological hori- 
zon the black hole would expand at the same rate as the 



universe (Zel'dovitch and Novikov 1967 Carr et al. 20101 



However, as Carr et al. (2010) explain, if the latter were true 



we would expect the existence of black holes that are sev- 
eral orders of magnitude larger than the black holes found in 
active galactic nuclei. Further studies of PBH lead to proof 
that they cannot expand at the same rate as the universe as 



shown in Carr and Hawkingj ([1974J) . 

More recent research continues the debate over the exis- 
tence of observable PBHs. Some authors show that the cal- 
culated accretion rates for simplified models are sufficiently 



high to transform a PBH to solar mass black holes ( Sivaram 



2001). As cited above, others make claims that when more 



complicated accretion models are considered, the mass ac- 
cretion rates for both stationary and moving PBHs are in- 
sufficient to appreciably modify the mass of the PBH except 



in very high Lorentz factor models (Custodio and Horvath 
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2002 Horvath and Custodio 2005). We show, using our own 
simplified model, that the accretion rates of moving PBH is 
insufficient for the PBH to reach solar masses, in agreement 



with the conclusions of Custodio and Horvath ( 2002 ) 



Traditional use of perfect radiation fluids is seen in cos- 
mological studies, where one would expect the equation of 
state to be ( |Miller||2005[ ), 



A 
2^ 



(1) 



where P is the pressure, p is the mass-energy density, and 
A is the cosmological constant. If we generalize this rela- 
tion between the fluid pressure, P, and energy density p, we 
havfQ 

^ = (r-i)(p-4Y (2) 



(3) 



2tt J 

This equation of state gives us a speed of sound, 

c 2 = ^=r-i 

dp 

Due to relativistic limits on the speed of sound we restrict 
our attention to V < 2. Furthermore, we are interested in 
treating the background matter as a fluid so we also restrict 
our attention to matter with F > 1. When F = 1 the back- 
ground fluid becomes a pressure-less dust. 

The equation of state |2| is valid in the very early uni- 
verse, ~ 10 _5 s after the Big Bang close to the time of the 
cosmological quark-hadron transition. This is deep inside 
the radiation-dominated era (.Miller 2005" I. During this time 
we expect the temperatures to be on the order of 10 12 K, 
mass-energy densities and horizon scales on the order of neu- 
tron star central densities, ~ 10 15 g/cm 3 , and radius ~ 10km 
respectively. The energy density of the early universe was 
dominated by relativistic particles and radiation. Thus equa- 
tion of state |2| is valid for any part of the radiation domi- 



nated phase, up to ~ 10 s after the Big Bang (Miller 2005) 



Consequently, during this time, we expect our equation of 
state to be valid, even on scales close to the black hole size 
~ 10 15 g. 

Our model assumes that the background fluid has uni- 
form density. The fluid itself is not sufficiently dense to ex- 
perience an appreciable self-gravity. We assume that for the 
duration of these simulations that the background space- 
time is stationary. This assumption is in agreement with 



results from Custodio and Horvath (20021, and allows us to 



treat the accretion problem using relativistic Bondi ( Michel 
1972) and relativistic Bondi-Hoyle ( Penner||2011 Font and 



Ibanez||1998a|b| |Font et al.||1999| |Petrich et al.||1989| |1988 1 

accretion models. Furthermore, this paper focuses on ax- 
isymmetric accretion onto an axisymmetric black hole, the 
next paper in this series studies a similar accretion system 
in an infinitely-thin disc model. 



In previous studies by Font and Ibanez (1998a I, an ul 



trarelativistic system was investigated and determined to 
reach a steady state. However, their physical model spec- 
ified the speed of sound in the background fluid to be 
c s = \/r — 1, while maintaining the usual ideal fluid en- 
thalpy per unit rest mass, 

h = 1 + e + P/p . (4) 



1 We used geometric units so that the speed of light is unity. 



Our model of the background fluid's equation of state en- 
forces the ultrarelativistic model such that the enthalpy den- 
sity has the form, 



h = p + P. 



(5) 



In our picture p does not include the rest mass density. 

For simplicity, our model assumes that the cosmolog- 
ical constant, A, is zero. While this constant is important 
for the dynamic study of gravitational collapse, in the early 
universe, A does not change the speed of sound Q and thus 
does not change the characteristic parameters in the system. 
Consequently, the cosmological constant is not expected to 
be dynamically important in our fixed spacetime model. 

To be consistent with "standard" radiation fluids we 
primarily investigate fluids with adiabatic index T = 4/3. 
However, we extend our model to investigate a range of adi- 
abatic constants and found similar behaviour in all cases. 
The parameters investigated may be found in table [T] For 
an ultrarelativistic fluid with T = 2 we refer the reader to 
|1988J, where a closed form solution is pre- 



Petrich et al. 



sented. We investigate three regions of parameter space, the 
subsonic regime < c s , the marginally supersonic regime 
«oo > c s , and the supersonic regime Voa > c s . As expected 
from previous studies the supersonic regimes contain the 
most interesting dynamics. 

Another interest of previous relativistic Bondi-Hoyle 



studies (Font and Ibanez 


1998a|b Font et al. 


1999 Pa- 


padopoulos and Font|1998 


Donmez et al.|2010l 1 


bcused on 



the search for an instability that would develop as the fluid 
accretes onto the black hole (Font and Ibanez 1998a). This 



instability is known as the "flip-flop" instability and is real- 
ized by the motion of the tail shock that rapidly moves from 
one side of the black hole to the other, eventually moving 
to the upstream side of the black hole destroying the stable 



accretion flow. Font and Ibanez (1998a) show that the evo- 



lution of the axisymmetric relativistic Bondi-Hoyle setup 
result in steady-state accretion rates for the parameters in- 
vestigated. Since any instabilities would be strongest on the 
event horizon, we calculate the accretion rates on the event 
horizon of the black hole. We further show, for the duration 
of our simulations, that instabilities were not found. 

To parameterise our fluid we use the ratio of the asymp- 
totic fluid speed to the asymptotic speed of sound, known 
as the asymptotic Mach number, 



(6) 



In analogy to the relativistic Lorentz factor W = 
1 — v J 2 , we define the sound speed Lorentz factor W c = 
J 2 (Font and Ibanez 



rections, the relativistic Mach 



1998a). Including these cor- 



number is defined as (Konigl 



1980), 



M 



R 



Wv 



M. 



(7) 



W c c B W c 

We further define the asymptotic relativistic Mach number 
as. 



M 



R 



Wg°ci 



Wo, 



Mc 



(8) 



The relativistic Mach number is the ratio of the proper speed 
of the fluid with the proper speed of sound in the fluid, 
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when this value is greater than one the flow is designated 
supersonic. When the Mach number is close to, but greater 
than one, we say that the flow is marginally supersonic. 

The rest of the paper will proceed as follows; in section[2] 
we present the coordinates used to study this problem. Sec- 
tion [3] discusses the equations of motion used in our ultra- 
relativistic system. Section [4] presents the initial conditions 
and boundary conditions used to perform the simulations. 
In section [5] we briefly cover the numerical methods devel- 
oped for the simulations. Section [6] is where we discuss the 
accretion profiles, while in sectionlTlwe discuss the accretion 
rates. In section [8] we present the flow morphology. Finally, 
section [9] contains our conclusions. 

In this paper we use geometric units, where G = c = 1 
with c the speed of light in vacuum, and G being New- 
ton's gravitational constant. We also use the convention that 
Greek indices are spacetime indices, while Latin indices are 
purely spatial indices. 



2 COORDINATES 



As in our earlier study ( Penner||2011 1, we are interested in 
the flow around a rotating black hole. The line element used 



is originally presented in Papadopoulos and Font (19981, 



ds 



2Mr \ 



df + 



4Mr 
A 



dtdr 



2Mr\ , 2 „ A 2Mr \ . - 
- I 1 + J dr 2 - 2a ( 1 + I sin 2 

. ,„ 2 4aMr sin 2 6 , , , 
+ Ad6 H dtd<j) 



+ (A- 



1 + 



2Mr 



Vs 



A J 



sin 6*) sin 



with 



A = r 2 + a 2 cos ( 



(9) 



(10) 



a is the dimensionless measure of the rotation rate of the 
black hole and is related to the angular momentum of the 
black hole via J = M 2 a, where M is the mass of the black 
hole. For the present study we set M = 1 without loss of 
generality. 



3 EQUATIONS OF MOTION 

We obtain the equations of motion for the ideal ultrarela- 
tivistic hydrodynamic (UHD) system by using the conser- 
vation of stress-energy. The stress-energy tensor used to de- 
scribe the UHD fluid is the same as that presented by |Neilsen| 
and Choptuik (20001. To close the system of equations, we 
use an equation of state ([2| with A = to relate the internal 
energy density to the fluid pressure. 



o 



(ii) 



where V M is the covariant derivative. In the UHD limit we 
write the hydrodynamic stress-energy tensor: 

T M " = hu^u + Pg^ (12) 

where h is the enthalpy density of the system defined in Q . 



We use the Anrowitt-Desner-Misner (ADM) 3 + 1 for- 
malism to re-express our system of partial differential equa- 



tions (York Jr. 19791, and use the Valencia formulation to 



describe our equations of motion. The equations of motion 



are presented in Neilsen and Choptuik (2000). We adapt 



their equations to a stationary spacetime background as seen 
below. 

The spatial momentum components are defined as, 



Sj = -n M 7„ J T'" / = hW 2 Vj 



and the energy, 



E = n, t n„T M " = hW 2 



(13) 



(14) 



where W is the relativistic Lorentz factor, Vj are the 3+1 
3-velocities, and P is the fluid pressure. 

This allows us to write the equations of motion as, 

d 



9 ^ d 




INITIALIZATION AND BOUNDARY 
CONDITIONS 



We use the method described by Font and Ibancz ( 19 98a|b I; 
Penner (2011 1 to initialize the hydrodynamic variables, with 



the exception of the baryon density po which is neglected in 
our model. 

The domain of integration for this study is defined by 



and < 6 ^ 7r. r n 



is determined in 



such a way that it will always fall inside the event horizon. 



The maximum radial domain, 



, was set to be sufficiently 



far from the event horizon that it would be effectively con- 
sidered infinity. If the radial domain is not set to be large 
enough, the bow shocks interfere with the outer, upstream, 
boundary conditions which quickly destroy the simulation. 

While previous studies assumed that a domain size 
r max = 50M was sufficiently large for ultrarelativistic sys- 
tems, we discovered that this is not true for our system. 
Using such a small radial domain resulted in a simulation 
that was marred with unphysical boundary effects. We in- 
vestigated this dependence and found that the formation of 
a steady state accretion flow was highly dependent on the 
radial distance between the black hole and domain edge. 
In particular we found that as the fluid becomes stiffer the 
larger the domain of integration needed to be. In this study 
we used regular grid spacing for both the radial and polar 
coordinates. In future studies we will investigate the use of 
a geometrically spaced grid. The geometric grid spacing is 
expected to provide higher accuracy near the event horizon 
of the black hole, while remaining numerically stable. 

The numerical treatment of the boundary conditions 
for this problem are the same as those presented in |Penner] 
( |2011[ ); |Font and Ibanez] ( |1998a[ ). 
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5 NUMERICAL METHODS 



The equations of motion (151 and (161 take on a general 
form, 

|V7Q+^^F*(Q) = ^S, (17) 

where Q are the conservative variables, F 1 denotes the flux, 
and S are the geometric source terms. To solve this sys- 
tem, we modified our high resolution shock capturing code 
used for the GRMHD (general relativistic magnetohydro- 



dynamic) study presented in Penner (20111. We changed 



the conservative variables, the primitive variable recovery 
scheme, and introduced a floor, as discussed below. Conver- 
gence tests may be found in App. [A] 

Parallelization was performed using the Parallel Adap- 
tive Mesh Refinement (PAMR). The PAMR infrastructure 



was developed by |Pretorius and Choptuik ( 2006 1 and was 



built on the message passing interface C software. For the 
numerical simulations presented for models U1-U4 we use 
a 400 x 160 grid. When we increased the domain size we 
adjusted the radial grid to maintain the same spatial res- 
olution. Simulations were performed using the Iridis clus- 
ter at the University of Southampton, the woodhen cluster 
at Princeton University, USA, and the WestGrid cluster of 
Canada. 



5.1 Primitive Variable Recovery 

Since we use a different set of primitive and conservative 
variables we describe the primitive variable recovery scheme. 
Our method used a modified version of Del Zanna's one pa- 
rameter inversion scheme ( Zanna et al.|2003 1 most recently 
investigated by |Noble et al.| ( |2006[ )" 
We define, 



Q = hW z 



(18) 



The first step of the primitive variable recovery is to 
solve the conservative equations for the term SI, as outlined 
below. 

We write the energy equation \14\ , 

(19) 



e = n - p, 



and using the momentum equation ( 13 1 we get 



|sf (n) = 7 1J SiSj = ar[i- 



1 

W 2 



(20) 



We solve (20) for W, 

w(n) = 



i - 



\sl 



(21) 



then re-formulate the equation of state, P — (T — l)p as 

fi(r- 1) 



P(fi) = 



vw 2 



(22) 



With the primitive variables defined as functions of Q. we use 
an iterative scheme to determine the values of the primitive 
variables. Using Eqn. ( |19[ ), we define /(fi), 

(23) 



/(H) =Q-P-E = 0. 



We then use Newton's method to solve the nonlinear equa- 

/(«') 



tion ( 23 1 for Q. 



i+i 



(24) 



(25) 



where I denotes the Newton iteration. To calculate df /dQ., 
we use the relations listed below, 

dCl 
dP 

dn 

dW 



= i-p' 

(w - 2Q.W) (r - 



rw 3 



(26) 
(27) 



= -W 



\s\ 2 



(28) 



dQ, " 

where the prime denotes a derivative with respect to fl. We 
perform this procedure at every time step, n. The Newton 
iteration must be initialized, so as an initial guess for the 

10 



value of fi we use the definition of Q, in Eqn. ( 18 \ and 
values for h, and W from the previous timestep, n — 1, 



n u = h n 



- 1 / / 1 n— 1 n— 1 \ 2 

/(l - JijVi «- 



, , (29) 
Given Q, we solve for the primitive variables. Using (1131) 



and ( 18 1 we express the velocity components as; 



(30) 



We then calculate the Lorentz factor W and pressure P from 



(21 1 and (22 l respectively. 



5.2 Floor 

There are instances where the code update routines produce 
conservative variables, where there are not any solutions for 
the primitive variables. This problem is independent of the 
inversion scheme used. In these instances, the primitive vari- 
ables are reconstructed using interpolated values as deter- 
mined by the primitive variables in surrounding cells. It is 
important to note that these instances were isolated to one 
or two cells in each time step, and also were only found in 
the downstream tail region of the evolution. Any effects due 
to the averaging were found to travel out of the domain, 
away from the black hole. It should also be noted that these 
numerical complications were restricted to the most extreme 
flow simulations. 



6 ACCRETION PROFILES 

In this section we describe the new and major features of 
each flow studied. This is done by presenting cross sections 
of the pressure accretion profiles found in Table [l] 

The first profiles we investigate are those for models 
U1-U4 with a — 0. We found, for the supersonic black holes 
with Doo = 0.6 and Voo — 0.9, that the pressure profiles along 
the axis of symmetry were dramatically different. In Fig. [l] 
we present the pressure profile for model U2 where we see 
a bow shock has formed. As we expect from previous hy- 
drodynamic studies ( Font and Ibanezjl 998a), the upstream 
pressure profile shows that the maximum upstream pressure 
is much smaller than the maximum pressure on the down- 
stream side of the black hole along the event horizon. The 
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Model 


r 


Vac 


^max 


* t hi 

lv± oo 


Ul 


4/3 


0.3 


100 


0.44 


U2 


- 


0.6 


- 


1.05 


U3 


- 


0.7 


- 


1.38 


U4 


- 


0.9 




2.91 


U5 


3/2 


0.6 




0.75 


U6 




0.7 




0.98 


U7 




0.71 




1.01 


U8 




0.9 




2.06 


U9 


1.1 


0.3 


: 


0.94 


U10 




0.32 




1.03 


Ull 




0.6 




2.27 


U12 


5/4 


0.4 


400 


0.86 


U13 




0.5 




1.02 


U14 




0.7 




3.62 



traditional tail shock seen in previous hydro dynamic studies 



Table 1. Table of parameters used for the axisymmetric systems 
studied in this paper. We have selected the velocities to capture 
the flows in the subsonic, marginally supersonic, and the super- 
sonic regime. We used the same velocity parameters for different 
values of the rotation parameter, a = 0.0, 0.5, 0.9, 0.99. 



unexpected feature of this flow is the presence of the bow 
shock, which was not obtained in previous studies due to the 
parameter range investigated. The detached or bow shock 
seen in the left plot in Fig. [I] appears to be a feature of the 
marginally supersonic flows in UHD systems. The detached 
shock persists when the black hole rotates, which may be 
seen in Fig. [4] 

When we increase the flow of the fluid past the black 
hole we find that the detached shock is prevented, as seen in 
Figs. [2] and [3] Model U3 begins to show signs that the tail 
shock will detach from the black hole; however, the shock 
stabilises near the front of the black hole, never making 
contact with the axis of symmetry. We see a snapshot of 
the stable tail shock in Fig. [9] 

When our simulation begins the fluid passes the black 
hole and focuses on the downstream side of the black hole, 
just as seen in previous studies (Font and Ibancz 1998a I. As 
the simulation proceeds the tail shock that forms widens, in- 
creasing the Mach cone opening angle. The fluid continues to 
flow past the black hole interacting with the pressure from 
fluid in the tail. The tail shock causes a back-pressure in 
the oncoming fluid. Depending on the value of the asymp- 
totic velocity one of two things happen. If the velocity is 
high enough, as seen in models U3 and U4 in the bottom 
panels of Fig. [5] the pressure in the tail shock will balance 
the fluid pressure from the oncoming fluid. Consequently, 
the tail shock stabilises. If, on the other hand, the asymp- 
totic fluid velocity is not sufficiently high as in models Ul 
and U2 in the top panels of Fig. [9] the pressure in the tail 
shock will continue to build, over-powering the oncoming 
fluid pressure. In this instance the tail shock will continue 
to widen, the angle of attachment will move to the upstream 
side of the black hole, make contact with the axis of sym- 
metry, and detach from the black hole. We see evidence of 
this when comparing models U2 and U3. While a bow shock 
does not form in U3, we see that the tail shock attaches to 
the black hole in the upstream region. The pressure from 
the oncoming fluid in model U3 is sufficient to balance the 
tail shock, and thus the tail shock persists. 

When we investigate black holes with more extreme 
asymptotic speeds such as seen in model U4 we recover the 



(Font and Ibafiez 1998a I. We present a snapshot of model 
U4 for all four spin parameters investigated in Fig. [9] It is 
clear that the tail shock persists. 

As we would expect for systems with tail shocks, when 
looking at the on-axis pressure profiles for models U3 and 
U4 in Figs. [2] and [3] respectively, the upstream pressure is 
dwarfed by the downstream pressure on the event horizon. 
By using horizon penetrating coordinates, we are sure this 
is not a coordinate effect as is cautioned by |Papadopoulos| 
and Font| ( [T998l ). 

When we investigated model U4 we obtain the expected 



tail shock presented in similar hydrodynamic studies (Font 
and Ibanez 1998a I 



We also investigated axisymmetric accretion onto a ro- 
tating black hole with the axis of rotation aligned with the 
axis of symmetry. The upstream pressure profiles were only 
marginally different for black holes with a = and a — 0.5 
as seen in the left plots in Fig. [5] When the spin rate 
was increased to a = 0.9 the upstream profile noticeably 
changed. For all spins investigated the upstream pressure 
remained low relative to the maximum downstream pres- 
sure. The downstream pressure profiles themselves were not 
appreciably altered, as is seen in the inset plot of Fig. [5] 
For a comparison of the effects of the black hole spin on 
the accretion profile we refer the reader to Fig. [10] where we 
present a snapshot of model U4 for all four spin parameters 
investigated. The spin rate of the black hole did not have a 
significant impact on the presence of a tail shock. The effect 
of the rotation of the black hole is more readily apparent in 
the angular cross section seen in Fig. [6] As the rotation rate 
of the black hole is increased the pressure on the axis of sym- 
metry decreases while the opening angle of the Mach cone 
increases. However, for the black hole spin rate to impact 
the type of shock exhibited the asymptotic fluid velocity 
will have to be slower than those investigated in this study. 
For the given parameters the black hole spin parameter will 
most likely effect the fluid flows with 0.6 < Voo < 0.7. We 
would also expect a parameter to exist that places the bow 
shock on the upstream side of the black hole. When we in- 
vestigate other adiabatic constants we see similar behaviour 
in all ultrarelativistic systems, as may be seen in [B] Mod- 
els of stiffer fluids have a faster speed of sound, and are 
expected to produce results with a wider Mach cone. Since 
our parameter survey searches for marginally supersonic pa- 
rameters as well as subsonic, we find that if we do not use 
a sufficiently large domain of integration that the tail shock 
will come in contact with the upstream outer boundary. This 
causes an unphysical pressure build-up at the outer domain 
that quickly dominates the fluid simulation. Thus for slower 
fluid flows we found that the domain of integration had to 
increase to capture the dynamics we are interested in. 



7 ACCRETION RATES 

For our ultrarelativistic system we measure the energy ac- 
cretion rate, 



E = 



gT 



sr 



> + / T» v T\udrd<t>, 



(31) 



where g is the determinant of the metric (j9J|. 
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Figure 1. Ultrarclativistic pressure profile upstream (top) and downstream (bottom) for model U2 with radial domain Teh ^ r ^ 100. 
To emphasize the bow shock in the upstream region we restrict our attention to the region —50 ^ r ^ —''eh- The pressure in the 
upstream region of the black hole is small by comparison to the pressure in the downstream region, or wake. The pressure profile in the 
upstream region indicates that there is a shock in the upstream region of the black hole. The above profiles were taken at t = 1000M. 
With the larger radial domain, we see that the pressure profile smoothly matches the upstream boundary conditions. The dashed line in 
the downstream profile emphasizes the maximum pressure of the fluid surrounding the black hole in the upstream side. 



Since we are simulating an ultrarelativistic fluid we can- 
not directly calculate the baryon mass accretion rate as per- 



formed in Font and Ibafiez ( |1998a|b[ ); Penner (20111. 

As in our previous study we measure the accretion rates 
at the Schwarzschild radius. We rescale the accretion rates 
with the equivalent Bondi accretion rates as determined 
from our simulations when the black hole is stationary rela- 
tive to the asymptotic background fluid. The energy accre- 
tion rates may be seen in Figs. [7] and [8] 

We see in figure [8] that the models with slower asymp- 
totic speeds (models U1,U2) reach a steady state around 
t — 200AZ and for the faster asymptotic speeds (models U3, 
U4) the steady state is reached near t = 600Af. The plots 
for the remaining models may be found in App.lB] We note, 
when comparing Fig. |B1| or Fig. |B3| with Fig. |B2| that for 
supersonic flows with a lower the adiabatic constant, the 
faster the system reaches the final steady state. 



8 FLOW MORPHOLOGY 

The results of the simulations using the parameters found 
in table [Tj are discussed here. We present the final state of 
model U1-U4 for both the spherically symmetric (a — 0) 
and axisymmetric (a ^ 0) black holes. Each of the param- 
eters studied here establish a steady state solution. As in 
the general hydrodynamic models ( Font and Ibanez|1998a l, 
we find that in supersonic flows, the upstream region of the 



flow is smooth, while in the downstream region there is the 
presence of a Mach cone. This is seen in the bottom right 
panel of Fig. [9] While for marginally supersonic flows such 
as seen in the bottom left panel of Fig. [9] the shock attaches 
to the black hole in the upstream region. In both models 
U3 and U4, the attached shock extends downstream, to the 
outer domain edge. 

The top two panels in Fig. [9] present the flow patterns 
during the steady state phase of the fluid flow for models Ul 
and U2 respectively. As expected for subsonic flows, seen in 
the top left panel, the pressure profile does not exhibit any 
shocks. The flow pattern changes when we increase the flow 
rate beyond the speed of sound. For r = 4/3, c s ~ 0.57 
so the asymptotic velocity Voo — 0.6 used for model U2 is 
marginally supersonic. The top right panel of Fig. [9] ex- 
hibits a detached/bow shock. Even in the presence of a bow 
shock the flow retains its steady state behaviour. Models U2 
and U3 appear to bound a special bow shock that persists 
along the event horizon. As seen in Fig. |10| when the black 
hole spins, with the spin axis aligned with the fluid axis of 
symmetry, the fluid profiles do not change appreciably rel- 
ative to a non-rotating (a — 0) black hole. This behaviour 
is expected since a detailed inspection of the pressure cross 
sections in Fig.[6]reveals that the cross-sections change only 
slightly as a function of spin. The upstream fluid experiences 
a slight increase in pressure while the downstream region 
remains similar for all spin parameters. Since these changes 
are only slight, we expect the overall fluid morphology to be 
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Figure 3. The pressure cross-section in the upstream region along the axis of symmetry plotted for model U4 in the region r^n ^ r ^ 50. 
As with model U3 in Fig.[2]the pressure is maximal on tbh- The downstream pressure on the event horizon dwarfs the upstream pressure 
on the event horizon by two orders of magnitude. 
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Figure 4. Ultrarelativistic pressure profiles for the rotating axisymmetric accretion problem for model U2, with a = 0.0, 0.5, 0.9, 0.99. 
We see that when the black hole rotates the bow shock persists. The upstream plot focuses on the region rj;H ^ ^ 20 to emphasize the 
bow shock region. We see the effects of rotation on the sharpness of the bow shock. As the spin of the black hole increases the resulting 
shock front begins to smooth. Since the shock front for a = 0.9 is only marginally sharper than that for a = 0.99, we suspect that the 
spin rate would need to exceed the maximal value allowed i.e. a > 1 to prevent the existence of a bow shock. Recall that if a > 1 wc 
have a naked singularity. 



similar, as seen in Fig. [6] The fact that the spin rate does 
not have a large impact on the fluid morphology is seen if 
we look back at the upstream pressure profiles for model U2 
in Fig. [4] where the bow shock persists with the spin turned 
on. We have seen that the shock opening angle is dependent 
on the black hole rotation rate, but for the parameters inves- 
tigated, we also find that the rotation rate does not impact 
the type of shock that forms. 



9 CONCLUSIONS 

We have shown that there is a relation between a system 
reaching a steady state and the size of the domain of inte- 
gration. The lower the adiabatic constant the larger the do- 
main needs to be prevent boundary effects from disrupting 
the simulation. In fact we also see that this is particularly 
important for subsonic flows which have detached shocks 
traveling upstream. 

Our simulations did not indicate the presence of any 
instabilities. This is likely due to the reduction in the com- 
plexity of the overall system. In the traditional black hole ac- 
cretion problems we are concerned with both the rest-mass 
interactions and the energy, while in our system the rest- 
mass contribution is negligible. Studies using axisymmetry 
tend to fail to capture instabilities due to the restrictive na- 
ture of the symmetry. This does not rule out the existence 
of instabilities similar to the standing accretion shock insta- 
bility (Fernandez 2010 1 or those presented by Donmez et al 



( 2010 1, where the instability timescale is very large, and con- 



sequently was not discovered using short-time simulations. 

It may be argued that running our simulations could be 
run for longer times, it becomes questionable how well we 
can approximate the accreting body as maintaining a fixed 
mass-energy, or how long before issues such as dynamical 
friction become important. Certainly, we expect that the 
assumption that the background fluid has constant pressure 
will become unreasonable over sufficiently long timescales. 

Our next study investigates the UHD system described 
above using an infinitely thin-disc accretion model, which re- 
laxes the axisymmetry. According to the studies by Donmez 



et al. (20101 using the usual hydrodynamic approximation, 



the thin-disc symmetry allows the formation of a flip-flop 
instability. 

The ultrarelativistic flow past a spherically symmet- 
ric black hole is in strong agreement with the flow past an 
axisymmetric/rotating black hole with the axis of rotation 
aligned with the direction of the fluid flow. The morphology 
of both geometric configurations is very similar and does 
not deviate too far from the results of previous studies of 
the full hydrodynamic evolution. Ultimately, most parame- 
ters surveyed resulted in a steady state solution. The only 
exception has been the r = 1.1 models which require a large 
amount of computational resources to reach a steady state. 
The accretion rates exhibited a dependence on the location 
of the computational boundary, such that a larger domain 
was essential to capture the proper flow patterns for some 
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Figure 5. Ultrarelativistic pressure profiles for the rotating axisymmctric accretion problem for model U4, with a = 0.0, 0.5, 0.9. We 
see that for black hole spin parameters a = and a = 0.5 the downstream (right) pressure profiles are only marginally distinguishable. 
As the spin is increased to a = 0.9, we see a much larger pressure in the upstream region of the black hole (left). The profiles on the 
downstream region become slightly sharper when the spin parameter is increased. 
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Figure 8. The energy accretion rates for models U2 (left) and U4 (right) for a = 0, 0.5, 0.9, 0.99. We see that the energy accreted by 
the black hole depends on the spin parameter a. As the spin parameter is increased the energy accretion drops significantly. The energy 
accretion rate is measured over a surface defined by the event horizon. As the spin parameter increases this surface decreases, thus we 
expect a decrease in the measured energy flux. 
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Figure 9. Pressure profile for the spherically symmetric evolution for models Ul with -M^> = 0.44 (top left), U2 with M.", = 1.05 (top 
right), U3 with = 1.38 (bottom left), and U4 with -MjJ, = 2.91 (bottom right), when a steady state is achieved. Model Ul was 

terminated at t = 400M due to boundary effects disrupting the simulation. We see that the black hole in model U4, traveling well above 
supersonic speeds, produces a readily apparent tail shock; however, when the flow is only marginally supersonic, as seen in model U2, 
we see the presence of a bow shock. If the flow is increased, as in model U3, we see that the tail shock persists, indicating there is some 
critical value where the bow shock remains fixed in the upstream region of the black hole along the event horizon. In the Newtonian 
studies, the presence of the bow shock was essential to the formation of what is known as the flip-flop instability. In our studies, the 
accretion rate on the event horizon remains constant and visualization of the flow remains constant, thus we conclude that the flow is 
steady despite the presence of a bow shock. In the first panel the flow is subsonic and so consequently, takes on a profile similar to that 
of Bondi accretion. 



adiabatic constants. We are currently investigating the ef- 
fects of magnetic fields on this same configuration. 

By comparing our accretion rates to that of the results 



Horvath and Custodio (2005 1, we also find that our mod- 



els of traveling black holes do not substantially increase their 
mass accretion rate. This leads us to agree with the con- 



clusions of Custodio and Horvath (2002) and Horvath and 



Custodio (2005), that a primordial black hole must sustain 



a very large boost to begin to accrete sufficient amounts of 
matter to completely avoid Hawking radiation. When com- 
paring fluid models for T — 1.1,5/4,4/3, and 3/2 in figures 
[BlJ[7||B3] and|B2j we find that the stiffer the fluid F — > 2, the 
faster the accretion rate. However, the Bondi-Hoyle accre- 
tion rates do not appear to increase beyond approximately 5 
to 7 times the Bondi accretion rate even for ultrarelativistic 
fluid flows with Voo > 0.9. We also find that the rotating 



black hole with spin parameter a — produces a maxi- 
mal accretion rate for all models, this is intuitive since the 
a — black hole has the largest radius, and consequently 
the largest "surface area". The fluids with adiabatic con- 
stant F — > 1 tend to accrete at a rate less than the Bondi 
accretion rate. 

To model more physically realistic accretion in the early 
universe future models will need to take into consideration 
the time dependent nature of the background energy den- 
sity, thus releasing the assumption of a uniform fluid back- 
ground. Furthermore, it will be beneficial for future stud- 
ies to consider general relativistic radiation hydrodynamics 



such as that presented by Park (2006) which may be used 
to approximate the radiative loses by the black hole due to 
Hawking radiation. 
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Figure 10. Pressure profile for the a = 0.0 (top left), a = 0.5 (top right), a = 0.9 (bottom left) and a = 0.99 (bottom 
right) axisymmetric evolution for model U4 at t = 400M. We see that a black hole traveling well above supersonic 
speeds produces a readily apparent tail shock. As the spin parameter of the black hole is increased we see that the 
point of contact for the tail with the event horizon of the black hole migrates upstream in agreement with the results 
presented in Fig. [6] All of these flows reach a steady state. Performing the same simulation using negative rotation 
rates produced identical pressure profiles and energy accretion rates. 
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APPENDIX A: CODE VERIFICATION 

Since the flow morphology is similar in all cases we present 
two convergence tests, we show the tests for Vao = 0.6 and 
Vao =0.9. These are seen in figures [AT] and |A2| The conver- 
gence values agree with those used by other authors. The dis- 
agreement between the true second order convergence value 
and my system comes from the special treatment of shock 
regions and local maxima where the system reduces to first 
order convergence. 
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Figure 6. Ultrarelativistic axisymmetric pressure profile 
P(1000, 4, 9) for a section through the tail shock at r = 4M for 
model U4 with spin parameters a = 0.0, 0.5, 0.9. The shock open- 
ing angle for black holes with spin greater than a = 0.5 is larger 
than that for smaller spin rates. 



APPENDIX B: ACCRETION RATES FOR 
DIFFERENT ADIABATIC CONSTANTS 

In this appendix we present the accretion rates for all the 
models studied in Table [T] We see that the softer fluid mod- 
els have a smaller accretion rate and that the subsonic mod- 
els for r = 1.1 have not reached a steady state, as seen in 
Fig. |B1| To capture the steady state solution for this model 
requires a much larger domain of integration, and longer 
time evolution. As the fluids stiffen the accretion rates in- 
crease, while the domain of integration may decrease. This 
trend is apparent when viewing figures |B"T| |B2[ [7] and |B3| 
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Figure Bl. The energy accretion rates for models U9 to Ull, T = 1.1, for all spin parameters. We see that the supersonic flow, v = 0.6, 
achieves a steady state solution early in the evolution; however, the subsonic flows have not reached a steady state. To capture the 
subsonic steady state solution we require the domain of integration to be extended much further than 400M, and we must run the 
simulations for longer times. 
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Figure B2. The energy accretion rates for models U5 to U8, T = 3/2. This model has a larger accretion rate, and reaches a steady state 
solution quicker than the softer fluid models. 
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Figure Al. A convergence test for model U2. We used a cross- 
section of the pressure variable in the upstream region (denoted 
by the negative radial direction) along the axis of symmetry. We 
see that the system is converging on a shock front upstream of the 
traveling black hole. Level 1 is 200 X 80, level 2 denotes 400 X 160, 
and level 3 denotes 800 X 320. In this convergence test we focused 
on the region r mln ^ r ^ 20 to emphasize the shock front. 



Figure A2. A convergence test for model U4. In this test we 
used a cross section along r = 4 around the black hole. We see 
that the solution is converging on a tail shock opening angle of 
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of approx. 107. Level 1 is 200 X 80, level 2 denotes 400 X 160, and 
level 3 denotes 800 X 320. 
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